Chapter 10 Does the antibiotic administration change the microbial community?

10.1 Alpha diversity

alpha_div_meta <- alpha_div %>%
  left_join(sample_metadata, by = join_by(sample == sample))%>%
  filter(time_point %in% c("Acclimation", "Antibiotics") ) 

label_map <- c(
  "richness" = "Species Richness",
  "neutral" = "Neutral Diversity",
  "phylogenetic" = "Phylogenetic Diversity")

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = "sample") %>%
  mutate(metric = factor(metric, levels = c("richness", "neutral", "phylogenetic")),
  time_point = factor(time_point, levels = c("Wild", "Acclimation"))) %>%
  filter(time_point %in% c("Wild", "Acclimation")) %>% 
  ggplot(aes(y = value, x = time_point, color=Population, fill=Population)) +
  geom_boxplot(width = 0.5,alpha=0.5, outlier.shape = NA, show.legend = FALSE) +
  geom_jitter(width = 0.2,show.legend = FALSE) +
  scale_color_manual(values=c('#008080', "#d57d2c")) +
  scale_fill_manual(values=c('#008080', "#d57d2c")) +
  facet_grid(metric ~ Population, scales = "free_y", labeller = labeller(metric = label_map, type = label_map))+
  theme_classic() +
  theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.text = element_text(size = 14),
    axis.title = element_text(size = 14),
    strip.text = element_text(size = 18, lineheight = 0.6),
      ) +
  ylab("Alpha diversity")

Richness

Modelq0GLMMNB <- glmer.nb(richness ~ time_point*Population+(1|individual), data = alpha_div_meta)
summary(Modelq0GLMMNB)

Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
 Family: Negative Binomial(3.3031)  ( log )
Formula: richness ~ time_point * Population + (1 | individual)
   Data: alpha_div_meta

     AIC      BIC   logLik deviance df.resid 
   453.0    464.6   -220.5    441.0       45 

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-1.45683 -0.56970  0.01616  0.42638  3.03485 

Random effects:
 Groups     Name        Variance Std.Dev.
 individual (Intercept) 0.1434   0.3787  
Number of obs: 51, groups:  individual, 27

Fixed effects:
                                     Estimate Std. Error z value Pr(>|z|)    
(Intercept)                            3.8038     0.1706  22.301  < 2e-16 ***
time_pointAntibiotics                 -1.2995     0.2143  -6.065 1.32e-09 ***
PopulationWarm                         0.6580     0.2833   2.323   0.0202 *  
time_pointAntibiotics:PopulationWarm  -0.1983     0.3580  -0.554   0.5797    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) tm_pnA PpltnW
tm_pntAntbt -0.539              
PopulatnWrm -0.603  0.324       
tm_pntAn:PW  0.335 -0.587 -0.553

Neutral

Model_neutral <- nlme::lme(fixed = neutral ~ time_point*Population, data = alpha_div_meta,
               random = ~ 1 | individual)
summary(Model_neutral)

Linear mixed-effects model fit by REML
  Data: alpha_div_meta 
       AIC      BIC    logLik
  357.0216 368.1225 -172.5108

Random effects:
 Formula: ~1 | individual
        (Intercept) Residual
StdDev:    3.334606 7.920952

Fixed effects:  neutral ~ time_point * Population 
                                         Value Std.Error DF   t-value p-value
(Intercept)                           19.30078  2.083103 25  9.265401   0e+00
time_pointAntibiotics                -12.24105  2.728975 22 -4.485585   2e-04
PopulationWarm                        25.30058  3.542048 25  7.142924   0e+00
time_pointAntibiotics:PopulationWarm -20.55961  4.732066 22 -4.344742   3e-04
 Correlation: 
                                     (Intr) tm_pnA PpltnW
time_pointAntibiotics                -0.655              
PopulationWarm                       -0.588  0.385       
time_pointAntibiotics:PopulationWarm  0.378 -0.577 -0.638

Standardized Within-Group Residuals:
       Min         Q1        Med         Q3        Max 
-1.8337660 -0.5250331 -0.1690052  0.5647389  2.1144650 

Number of Observations: 51
Number of Groups: 27 

Phylogenetic

Model_phylo <- nlme::lme(fixed = phylogenetic ~ time_point*Population, data = alpha_div_meta,
               random = ~ 1 | individual)
summary(Model_phylo)

Linear mixed-effects model fit by REML
  Data: alpha_div_meta 
       AIC     BIC    logLik
  214.9941 226.095 -101.4971

Random effects:
 Formula: ~1 | individual
        (Intercept) Residual
StdDev:   0.8604909  1.70013

Fixed effects:  phylogenetic ~ time_point * Population 
                                         Value Std.Error DF   t-value p-value
(Intercept)                           5.426628 0.4616166 25 11.755705  0.0000
time_pointAntibiotics                -1.101856 0.5866694 22 -1.878155  0.0737
PopulationWarm                        1.088797 0.7851891 25  1.386669  0.1778
time_pointAntibiotics:PopulationWarm -0.375271 1.0172677 22 -0.368901  0.7157
 Correlation: 
                                     (Intr) tm_pnA PpltnW
time_pointAntibiotics                -0.635              
PopulationWarm                       -0.588  0.374       
time_pointAntibiotics:PopulationWarm  0.366 -0.577 -0.618

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-1.79893698 -0.50285915 -0.07891184  0.49006037  1.85177334 

Number of Observations: 51
Number of Groups: 27 

10.2 Beta diversity

samples_to_keep <- sample_metadata %>%
  filter(time_point %in% c("Acclimation", "Antibiotics")) %>% 
  select(Tube_code) %>% 
  pull()
subset_meta <- sample_metadata %>%
  filter(time_point %in% c("Acclimation", "Antibiotics"))

Richness

richness <- as.matrix(beta_q0n$S)
richness <- as.dist(richness[rownames(richness) %in% samples_to_keep,
               colnames(richness) %in% samples_to_keep])

betadisper(richness, subset_meta$time_point) %>% permutest(.) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq  Mean Sq      F N.Perm Pr(>F)    
Groups     1 0.046108 0.046108 12.545    999  0.001 ***
Residuals 49 0.180100 0.003676                         
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

adonis2(richness ~ time_point*Population,
        data = subset_meta %>% arrange(match(Tube_code,labels(richness))),
        permutations = 999,
        strata = subset_meta %>% arrange(match(Tube_code,labels(richness))) %>% pull(individual),
        by="terms") %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
time_point	1	2.0239655	0.09924051	6.197237	0.001
Population	1	2.1298900	0.10443427	6.521570	0.001
time_point:Population	1	0.8908892	0.04368271	2.727839	0.004
Residual	47	15.3498055	0.75264251	NA	NA
Total	50	20.3945503	1.00000000	NA	NA

Neutral

neutral <- as.matrix(beta_q1n$S)
neutral <- as.dist(neutral[rownames(neutral) %in% samples_to_keep,
               colnames(neutral) %in% samples_to_keep])

betadisper(neutral, subset_meta$time_point) %>% permutest(.) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)   
Groups     1 0.06612 0.066120 12.357    999  0.002 **
Residuals 49 0.26219 0.005351                        
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

adonis2(neutral ~ time_point*Population,
        data = subset_meta %>% arrange(match(Tube_code,labels(neutral))),
        permutations = 999,
        strata = subset_meta %>% arrange(match(Tube_code,labels(neutral))) %>% pull(individual),
        by="terms") %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
time_point	1	2.1502445	0.1099662	7.412706	0.001
Population	1	2.8814616	0.1473616	9.933488	0.001
time_point:Population	1	0.8884255	0.0454352	3.062739	0.003
Residual	47	13.6335491	0.6972370	NA	NA
Total	50	19.5536807	1.0000000	NA	NA

Phylogenetic

phylo <- as.matrix(beta_q1p$S)
phylo <- as.dist(phylo[rownames(phylo) %in% samples_to_keep,
               colnames(phylo) %in% samples_to_keep])
betadisper(phylo, subset_meta$time_point) %>% permutest(.)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq Mean Sq     F N.Perm Pr(>F)    
Groups     1 0.82840 0.82840 60.21    999  0.001 ***
Residuals 49 0.67417 0.01376                        
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

adonis2(phylo ~ time_point*Population,
        data = subset_meta %>% arrange(match(Tube_code,labels(phylo))),
        permutations = 999,
        strata = subset_meta %>% arrange(match(Tube_code,labels(phylo))) %>% pull(individual),
        by="terms") %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
time_point	1	2.0945462	0.22473442	16.243398	0.001
Population	1	0.7920735	0.08498556	6.142603	0.001
time_point:Population	1	0.3729411	0.04001473	2.892192	0.031
Residual	47	6.0605345	0.65026529	NA	NA
Total	50	9.3200953	1.00000000	NA	NA

dbRDA

#Richness
cca_ord <- capscale(formula = richness ~ subset_meta$time_point* subset_meta$Population)
CAP_df <- as.data.frame(vegan::scores(cca_ord, display = "sites")) %>%
  rownames_to_column('Tube_code') %>%
  left_join(subset_meta, by = 'Tube_code') %>%
  column_to_rownames('Tube_code')%>%
  mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(CAP2, na.rm = TRUE))

biplot_scores <- as.data.frame(vegan::scores(cca_ord, display = "bp")) %>%
  rownames_to_column("Variable")
biplot_scores$Variable <- recode(
  biplot_scores$Variable,
  "subset_meta$time_pointAntibiotics" = "Antibiotics",
  "subset_meta$PopulationWarm" = "Warm population",
  "subset_meta$time_pointAntibiotics:subset_meta$PopulationWarm" =
    "Interaction Wam population"
)

CAP_df %>%
  group_by(Population, time_point) %>%
  mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(CAP2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(., aes(x=CAP1,y=CAP2, color=Population,shape = time_point)) +
  scale_color_manual(values=c('#008080', "#d57d2c")) +
  scale_fill_manual(values=c('#00808050', "#d57d2c50")) +
  geom_point(size=2) +
  geom_hline(yintercept = 0, linetype = "dotted") + 
  geom_vline(xintercept = 0, linetype = "dotted") +
  geom_segment(aes(x=x_cen, y=y_cen, xend=CAP1, yend=CAP2), alpha=0.2) +
  geom_segment(data = biplot_scores, aes(x = 0, y = 0, xend = CAP1, yend = CAP2),
               inherit.aes = FALSE, arrow = arrow(length = unit(0.2, "cm")), color = "black") +
  geom_text(data = biplot_scores, aes(x = CAP1, y = CAP2, label = Variable),
            inherit.aes = FALSE, color = "black", vjust = -0.5, hjust = 0.5)+
  theme_classic()

#Neutral
cca_ord <- capscale(formula = neutral ~ subset_meta$time_point* subset_meta$Population)
CAP_df <- as.data.frame(vegan::scores(cca_ord, display = "sites")) %>%
  rownames_to_column('Tube_code') %>%
  left_join(subset_meta, by = 'Tube_code') %>%
  column_to_rownames('Tube_code')%>%
  mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(CAP2, na.rm = TRUE))

biplot_scores <- as.data.frame(vegan::scores(cca_ord, display = "bp")) %>%
  rownames_to_column("Variable")
biplot_scores$Variable <- recode(
  biplot_scores$Variable,
  "subset_meta$time_pointAntibiotics" = "Antibiotics",
  "subset_meta$PopulationWarm" = "Warm population",
  "subset_meta$time_pointAntibiotics:subset_meta$PopulationWarm" =
    "Interaction Wam population"
)

CAP_df %>%
  group_by(Population, time_point) %>%
  mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(CAP2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(., aes(x=CAP1,y=CAP2, color=Population,shape = time_point)) +
  scale_color_manual(values=c('#008080', "#d57d2c")) +
    scale_fill_manual(values=c('#00808050', "#d57d2c50"))+
  geom_point(size=2) +
  geom_hline(yintercept = 0, linetype = "dotted") + 
  geom_vline(xintercept = 0, linetype = "dotted") +
  geom_segment(aes(x=x_cen, y=y_cen, xend=CAP1, yend=CAP2), alpha=0.2) +
  geom_segment(data = biplot_scores, aes(x = 0, y = 0, xend = CAP1, yend = CAP2),
               inherit.aes = FALSE, arrow = arrow(length = unit(0.2, "cm")), color = "black") +
  geom_text(data = biplot_scores, aes(x = CAP1, y = CAP2, label = Variable),
            inherit.aes = FALSE, color = "black", vjust = -0.5, hjust = 0.5)+
  theme_classic()

#Phylogenetic
cca_ord <- capscale(formula = phylo ~ subset_meta$time_point* subset_meta$Population)
CAP_df <- as.data.frame(vegan::scores(cca_ord, display = "sites")) %>%
  rownames_to_column('Tube_code') %>%
  left_join(subset_meta, by = 'Tube_code') %>%
  column_to_rownames('Tube_code')%>%
  mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(CAP2, na.rm = TRUE))

biplot_scores <- as.data.frame(vegan::scores(cca_ord, display = "bp")) %>%
  rownames_to_column("Variable")
biplot_scores$Variable <- recode(
  biplot_scores$Variable,
  "subset_meta$time_pointAntibiotics" = "Antibiotics",
  "subset_meta$PopulationWarm" = "Warm population",
  "subset_meta$time_pointAntibiotics:subset_meta$PopulationWarm" =
    "Interaction Wam population"
)

CAP_df %>%
  group_by(type, time_point) %>%
  mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(CAP2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(., aes(x=CAP1,y=CAP2, color=Population,shape = time_point)) +
  scale_color_manual(values=c('#008080', "#d57d2c")) +
  scale_fill_manual(values=c('#00808050', "#d57d2c50")) +
  geom_point(size=2) +
  geom_hline(yintercept = 0, linetype = "dotted") + 
  geom_vline(xintercept = 0, linetype = "dotted") +
  geom_segment(aes(x=x_cen, y=y_cen, xend=CAP1, yend=CAP2), alpha=0.2) +
  geom_segment(data = biplot_scores, aes(x = 0, y = 0, xend = CAP1, yend = CAP2),
               inherit.aes = FALSE, arrow = arrow(length = unit(0.2, "cm")), color = "black") +
  geom_text(data = biplot_scores, aes(x = CAP1, y = CAP2, label = Variable),
            inherit.aes = FALSE, color = "black", vjust = -0.5, hjust = 0.5)+
  theme_classic()

10.3 Differences in MAGs and functional capacity

10.3.1 Cold population

10.3.1.1 Structural zeros

acclimation_samples <- sample_metadata %>%
  filter(Population == "Cold" & time_point == "Acclimation") %>% 
  dplyr::select(sample) %>%
  pull()

antibiotics_samples <- sample_metadata %>%
  filter(Population == "Cold" & time_point == "Antibiotics") %>% 
  dplyr::select(sample) %>% pull()

subset_meta <- sample_metadata %>%
  filter(Population == "Cold" & time_point %in% c("Antibiotics","Acclimation"))

structural_zeros <- genome_counts_filt %>% 
  select(one_of(c("genome",subset_meta$sample))) %>%
  filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
  select(genome, where(~ is.numeric(.) && sum(.) > 0)) %>% 
   rowwise() %>% 
   mutate(all_zeros_acclimation = all(c_across(all_of(acclimation_samples)) == 0)) %>% 
   mutate(all_zeros_antibiotics = all(c_across(all_of(antibiotics_samples)) == 0)) %>% 
   mutate(average_acclimation = mean(c_across(all_of(acclimation_samples)), na.rm = TRUE)) %>% 
   mutate(average_antibiotics = mean(c_across(all_of(antibiotics_samples)), na.rm = TRUE)) %>% 
   filter(all_zeros_acclimation == TRUE || all_zeros_antibiotics==TRUE)  %>% 
   mutate(present = case_when(
      all_zeros_acclimation & !all_zeros_antibiotics ~ "antibiotics",
      !all_zeros_acclimation & all_zeros_antibiotics ~ "acclimation",
      !all_zeros_acclimation & !all_zeros_antibiotics ~ "None",
      TRUE ~ NA_character_
    )) %>%
   mutate(average = ifelse(present == "acclimation", average_acclimation, average_antibiotics)) %>%
   dplyr::select(genome, present, average) %>%
   left_join(genome_metadata, by=join_by(genome==genome)) %>%
   arrange(present,-average)

Acclimation

structural_zeros %>% 
  filter(present=="acclimation") %>% 
  count(phylum, name = "sample_count") %>%
  arrange(desc(sample_count)) 

# A tibble: 9 × 2
# Rowwise: 
  phylum               sample_count
  <chr>                       <int>
1 p__Bacillota_A                 43
2 p__Bacteroidota                15
3 p__Bacillota                    8
4 p__Pseudomonadota               8
5 p__Cyanobacteriota              3
6 p__Verrucomicrobiota            2
7 p__Bacillota_B                  1
8 p__Bacillota_C                  1
9 p__Fusobacteriota               1

structural_zeros %>% 
  filter(present=="acclimation") %>% 
  select(1,5:10)

# A tibble: 82 × 7
# Rowwise: 
   genome                phylum               class               order                 family                 genus          species
   <chr>                 <chr>                <chr>               <chr>                 <chr>                  <chr>          <chr>  
 1 AH1_2nd_1:bin_000006  p__Bacillota_A       c__Clostridia       o__Lachnospirales     f__Lachnospiraceae     g__Kineothrix  s__    
 2 AH1_2nd_19:bin_000005 p__Bacillota_A       c__Clostridia       o__Lachnospirales     f__Lachnospiraceae     g__Roseburia   s__    
 3 AH1_2nd_11:bin_000008 p__Verrucomicrobiota c__Verrucomicrobiae o__Verrucomicrobiales f__Akkermansiaceae     g__Akkermansia s__    
 4 AH1_2nd_2:bin_000003  p__Bacillota_A       c__Clostridia       o__Lachnospirales     f__Lachnospiraceae     g__JAAYNV01    s__    
 5 AH1_2nd_10:bin_000007 p__Bacillota_A       c__Clostridia       o__Lachnospirales     f__                    g__            s__    
 6 AH1_2nd_15:bin_000043 p__Bacteroidota      c__Bacteroidia      o__Bacteroidales      f__Bacteroidaceae      g__Bacteroides s__    
 7 AH1_2nd_10:bin_000041 p__Bacteroidota      c__Bacteroidia      o__Bacteroidales      f__Rikenellaceae       g__Alistipes   s__    
 8 AH1_2nd_8:bin_000007  p__Bacillota         c__Bacilli          o__Erysipelotrichales f__Erysipelotrichaceae g__Dielma      s__    
 9 AH1_2nd_14:bin_000015 p__Bacillota_A       c__Clostridia       o__Christensenellales f__UBA3700             g__            s__    
10 AH1_2nd_5:bin_000001  p__Bacillota_C       c__Negativicutes    o__Selenomonadales    f__                    g__            s__    
# ℹ 72 more rows

Antibiotics

structural_zeros %>% 
  filter(present=="antibiotics") %>% 
  select(1,5:10)

# A tibble: 4 × 7
# Rowwise: 
  genome                phylum            class                  order                 family                genus         species            
  <chr>                 <chr>             <chr>                  <chr>                 <chr>                 <chr>         <chr>              
1 AH1_2nd_7:bin_000003  p__Pseudomonadota c__Gammaproteobacteria o__Enterobacterales   f__Enterobacteriaceae g__Proteus    s__Proteus cibarius
2 LI1_2nd_7:bin_000001  p__Bacillota_A    c__Clostridia          o__Clostridiales      f__Clostridiaceae     g__Sarcina    s__                
3 AH1_2nd_7:bin_000055  p__Bacillota      c__Bacilli             o__Mycoplasmatales    f__Mycoplasmoidaceae  g__Ureaplasma s__                
4 AH1_2nd_13:bin_000025 p__Bacillota_A    c__Clostridia          o__Christensenellales f__UBA3700            g__           s__                

Community elements differences in structural zeros

element_gift <- GIFTs_elements_community %>% 
  as.data.frame() %>% 
  rownames_to_column(., "sample") %>% 
  left_join(., sample_metadata[c(8,9)], by=join_by("sample"=="sample"))

uniqueGIFT_db<- unique(GIFT_db[c(2,4,5,6)]) %>% unite("Function",Function:Element, sep= "_", remove=FALSE)

element_gift %>%
    pivot_longer(-c(sample,time_point), names_to = "trait", values_to = "value") %>%
    group_by(trait) %>%
    summarise(p_value = wilcox.test(value ~ time_point)$p.value) %>%
    mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
    filter(p_adjust < 0.05)

# A tibble: 0 × 3
# ℹ 3 variables: trait <chr>, p_value <dbl>, p_adjust <dbl>

10.3.1.2 Differential abundance analysis: composition

phylo_samples <- sample_metadata %>%
  filter(Population == "Cold" & time_point %in% c("Acclimation", "Antibiotics") )%>% 
  column_to_rownames("sample") %>% 
  sample_data()
phylo_genome <- genome_counts_filt %>% 
  filter(!genome %in% structural_zeros$genome) %>% 
  select(one_of(c("genome",rownames(phylo_samples)))) %>%
  filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
  select(genome, where(~ is.numeric(.) && sum(.) > 0)) %>%
  column_to_rownames("genome") %>% 
  otu_table(., taxa_are_rows = TRUE)
phylo_taxonomy <- genome_metadata %>% 
  filter(genome %in% rownames(phylo_genome)) %>% 
  column_to_rownames("genome") %>% 
  dplyr::select(domain,phylum,class,order,family,genus,species) %>% 
  as.matrix() %>% 
  tax_table()

physeq_genome_filtered <- phyloseq(phylo_genome, phylo_taxonomy, phylo_samples)

ancom_rand_output_antib_accl_2803 = ancombc2(data = physeq_genome_filtered, 
                  assay_name = "counts", 
                  tax_level = NULL,
                  fix_formula = "time_point",
                  rand_formula = "(1|individual)",
                  p_adj_method = "holm", 
                  pseudo_sens = TRUE,
                  prv_cut =0, 
                  lib_cut = 0, 
                  s0_perc = 0.05,
                  group = NULL, 
                  struc_zero = FALSE, 
                  neg_lb = FALSE,
                  alpha = 0.05, 
                  n_cl = 2, 
                  verbose = TRUE,
                  global = FALSE, 
                  pairwise = FALSE, 
                  dunnet = FALSE, 
                  trend = FALSE,
                  iter_control = list(tol = 1e-5, max_iter = 20, verbose = FALSE),
                  em_control = list(tol = 1e-5, max_iter = 100),
                  lme_control = lme4::lmerControl(),
                  mdfdr_control = list(fwer_ctrl_method = "holm", B = 100), 
                  trend_control = NULL)

save(ancom_rand_output_antib_accl_2803, 
     file = "data/ancombc_antib_accl_2803.Rdata")

load("data/ancombc_antib_accl_2803.Rdata")

taxonomy <- data.frame(physeq_genome_filtered@tax_table) %>%
  rownames_to_column(., "taxon") %>%
  mutate_at(vars(order, phylum, family, genus, species), ~ str_replace(., "[dpcofgs]__", ""))

ancombc_rand_table_mag <- ancom_rand_output_antib_accl_2803$res %>%
  dplyr::select(taxon, lfc_time_pointAntibiotics, p_time_pointAntibiotics) %>%
  filter(p_time_pointAntibiotics < 0.05) %>%
  dplyr::arrange(p_time_pointAntibiotics) %>%
  merge(., taxonomy, by="taxon") %>%
  mutate_at(vars(phylum, species), ~ str_replace(., "[dpcofgs]__", ""))%>%
  dplyr::arrange(lfc_time_pointAntibiotics)

colors_alphabetic <- read_tsv("https://raw.githubusercontent.com/earthhologenome/EHI_taxonomy_colour/main/ehi_phylum_colors.tsv") %>%
  mutate_at(vars(phylum), ~ str_replace(., "[dpcofgs]__", ""))  %>%
  right_join(taxonomy, by=join_by(phylum == phylum)) %>%
  dplyr::select(phylum, colors) %>%
  mutate(colors = str_c(colors, "80"))  %>% #add 80% alpha
    unique() %>%
    dplyr::arrange(phylum)

tax_table <- as.data.frame(unique(ancombc_rand_table_mag$phylum))
  
colnames(tax_table)[1] <- "phylum"
tax_color <- merge(tax_table, colors_alphabetic, by="phylum")%>%
    dplyr::arrange(phylum) %>%
    dplyr::select(colors) %>%
    pull()

ancombc_rand_table_mag%>%
      mutate(genome=factor(taxon,levels=ancombc_rand_table_mag$taxon)) %>%
ggplot(., aes(x=lfc_time_pointAntibiotics, y=forcats::fct_reorder(genome,lfc_time_pointAntibiotics), fill=phylum)) + 
  geom_col() + 
  scale_fill_manual(values=tax_color) + 
  geom_hline(yintercept=0) + 
  theme(panel.background = element_blank(),
        axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
        axis.text.x = element_text(size = 12),
        axis.text.y = element_text(size = 8),
        axis.title = element_text(size = 14, face = "bold"),
        legend.text = element_text(size = 12),
        legend.title = element_text(size = 14, face = "bold"),
        legend.position = "right", legend.box = "vertical")+
  xlab("log2FoldChange") + 
  ylab("Species")+
  guides(fill=guide_legend(title="Phylum"))

MAGs enriched in acclimation on the left side and in antibiotics on the right side

Phyla of the significant MAGs

Acclimation

ancombc_rand_table_mag%>%
  filter(lfc_time_pointAntibiotics<0)  %>% 
  count(phylum, name = "sample_count") %>%
  arrange(desc(sample_count))  

            phylum sample_count
1     Bacteroidota           15
2      Bacillota_A            4
3        Bacillota            1
4 Campylobacterota            1

ancombc_rand_table_mag%>%
  filter(lfc_time_pointAntibiotics<0)  %>% 
  select(phylum, genus, species)

             phylum           genus species
1  Campylobacterota  Helicobacter_J        
2       Bacillota_A                        
3      Bacteroidota     Bacteroides        
4      Bacteroidota     Odoribacter        
5         Bacillota                        
6      Bacteroidota     Bacteroides        
7      Bacteroidota     Bacteroides        
8      Bacteroidota     Bacteroides        
9      Bacteroidota     Phocaeicola        
10     Bacteroidota     Odoribacter        
11     Bacteroidota     Phocaeicola        
12     Bacteroidota     Bacteroides        
13     Bacteroidota     Bacteroides        
14     Bacteroidota Parabacteroides        
15     Bacteroidota                        
16      Bacillota_A                        
17     Bacteroidota       Alistipes        
18      Bacillota_A                        
19      Bacillota_A                        
20     Bacteroidota       Alistipes        
21     Bacteroidota     Odoribacter        

Antibiotics

ancombc_rand_table_mag%>%
  filter(lfc_time_pointAntibiotics>0)  %>% 
  count(phylum, name = "sample_count") %>%
  arrange(desc(sample_count))  

             phylum sample_count
1       Bacillota_A            2
2    Pseudomonadota            2
3 Verrucomicrobiota            1

ancombc_rand_table_mag%>%
  filter(lfc_time_pointAntibiotics>0)  %>% 
  select(phylum, genus, species)

             phylum          genus species
1       Bacillota_A     MGBC116941        
2    Pseudomonadota                       
3       Bacillota_A Intestinimonas        
4 Verrucomicrobiota                       
5    Pseudomonadota                       

10.3.1.3 Differences in the functional capacity

sample_sub <- sample_metadata %>%
  filter(Population == "Cold" & time_point %in% c("Acclimation", "Antibiotics"))
genome_counts_filtered <- genome_counts_filt %>% 
    select(one_of(c("genome",sample_sub$sample))) 
  
GIFTs_elements <- to.elements(genome_gifts,GIFT_db)
GIFTs_elements_filtered <- GIFTs_elements[rownames(GIFTs_elements) %in% genome_counts_filtered$genome,]
GIFTs_elements_filtered <- as.data.frame(GIFTs_elements_filtered) %>% 
  select_if(~ !is.numeric(.) || sum(.) != 0)

GIFTs_functions <- to.functions(GIFTs_elements_filtered,GIFT_db)

GIFTs_domains <- to.domains(GIFTs_functions,GIFT_db)

genome_counts_filtered_row <- genome_counts_filtered %>%
  mutate_at(vars(-genome),~./sum(.)) %>% 
  column_to_rownames(., "genome") 

GIFTs_elements_community <- to.community(GIFTs_elements_filtered,genome_counts_filtered_row,GIFT_db)
GIFTs_functions_community <- to.community(GIFTs_functions,genome_counts_filtered_row,GIFT_db)
GIFTs_domains_community <- to.community(GIFTs_domains,genome_counts_filtered_row,GIFT_db)

10.3.1.3.1 Community elements differences:

MCI

GIFTs_elements_community %>%
  rowMeans() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = join_by(sample == sample)) %>%
  group_by(time_point) %>%
  summarise(MCI = mean(value), sd = sd(value))

# A tibble: 2 × 3
  time_point    MCI     sd
  <chr>       <dbl>  <dbl>
1 Acclimation 0.331 0.0319
2 Antibiotics 0.244 0.132 

MCI <- GIFTs_elements_community %>%
  rowMeans() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = join_by(sample == sample)) 

shapiro.test(MCI$value)

    Shapiro-Wilk normality test

data:  MCI$value
W = 0.93705, p-value = 0.05038

wilcox.test(value ~ time_point, data=MCI)

    Wilcoxon rank sum exact test

data:  value by time_point
W = 209, p-value = 0.02612
alternative hypothesis: true location shift is not equal to 0

element_gift <- GIFTs_elements_community %>% 
  as.data.frame() %>% 
  rownames_to_column(., "sample") %>% 
  left_join(., sample_metadata[c(8,9)], by=join_by("sample"=="sample"))

uniqueGIFT_db<- unique(GIFT_db[c(2,4,5,6)]) %>% unite("Function",Function:Element, sep= "_", remove=FALSE)

significant_elements <- element_gift %>%
    pivot_longer(-c(sample,time_point), names_to = "trait", values_to = "value") %>%
    group_by(trait) %>%
    summarise(p_value = wilcox.test(value ~ time_point)$p.value) %>%
    mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
    filter(p_adjust < 0.05)%>%
  left_join(.,uniqueGIFT_db[c(1,3)],by = join_by(trait == Code_element))

element_gift_t <- element_gift  %>% 
  t() %>%
  row_to_names(row_number = 1) %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric)  %>%
  rownames_to_column(., "trait")

element_gift_filt <- subset(element_gift_t, trait %in% significant_elements$trait) %>% 
  t() %>%
  row_to_names(row_number = 1) %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric)  %>%
  rownames_to_column(., "sample")%>% 
  left_join(., sample_metadata[c(8,9)], by = join_by(sample == sample))

difference_table <- element_gift_filt %>%
  select(-sample) %>%
  group_by(time_point) %>%
  summarise(across(everything(), mean)) %>%
  t() %>%
  row_to_names(row_number = 1) %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric) %>%
  rownames_to_column(., "Elements") %>%
  left_join(.,uniqueGIFT_db[c(1,3,4)],by = join_by(Elements == Code_element)) %>% 
  arrange(Function) %>% 
  mutate(Difference=Acclimation-Antibiotics)%>% 
  mutate(group_color = ifelse(Difference <0, "Antibiotics","Acclimation")) 

difference_table %>%
  ggplot(aes(x=forcats::fct_reorder(Function,Difference), y=Difference, fill=group_color)) + 
  geom_col() +
  scale_fill_manual(values=c('#008080',"#003a45")) +
  geom_hline(yintercept=0) + 
  coord_flip()+
  theme(axis.text = element_text(size = 10),
        axis.title = element_text(size = 12),
        legend.position = "right", 
        panel.background = element_blank(),
          panel.grid.major = element_line(size = 0.15, linetype = 'solid',
                                colour = "grey"))+
  xlab("Microbial Functional Capacity") + 
  ylab("Mean difference")+
  labs(fill="Time_point")

10.3.1.3.2 Community functions differences

MCI

GIFTs_functions_community %>%
  rowMeans() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = join_by(sample == sample)) %>%
  group_by(time_point) %>%
  summarise(MCI = mean(value), sd = sd(value))

# A tibble: 2 × 3
  time_point    MCI     sd
  <chr>       <dbl>  <dbl>
1 Acclimation 0.329 0.0319
2 Antibiotics 0.254 0.123 

MCI <- GIFTs_functions_community %>%
  rowMeans() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = join_by(sample == sample)) 

shapiro.test(MCI$value)

    Shapiro-Wilk normality test

data:  MCI$value
W = 0.95036, p-value = 0.1259

wilcox.test(value ~ time_point, data=MCI)

    Wilcoxon rank sum exact test

data:  value by time_point
W = 213, p-value = 0.01774
alternative hypothesis: true location shift is not equal to 0

function_gift <- GIFTs_functions_community %>% 
  as.data.frame() %>% 
  rownames_to_column(., "sample") %>% 
  merge(., sample_metadata[c(8,9)], by="sample")

unique_funct_db<- GIFT_db[c(3,4,5)] %>% 
  distinct(Code_function, .keep_all = TRUE)

significant_functional <- function_gift %>%
  pivot_longer(-c(sample,time_point), names_to = "trait", values_to = "value") %>%
  group_by(trait) %>%
  summarise(p_value = wilcox.test(value ~ time_point)$p.value) %>%
  mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
  filter(p_adjust < 0.05)%>%
  left_join(.,unique_funct_db[c(1,3)],by = join_by(trait == Code_function))

Code_function	Acclimation	Antibiotics	Function	Difference	time_point
B06	0.68059020	0.47977960	Organic anion biosynthesis	0.20081060	Acclimation
D02	0.38530780	0.19353530	Polysaccharide degradation	0.19177250	Acclimation
S03	0.27145162	0.08296879	Spore	0.18848283	Acclimation
B02	0.59930820	0.42536320	Amino acid biosynthesis	0.17394500	Acclimation
B01	0.84159250	0.67849220	Nucleic acid biosynthesis	0.16310030	Acclimation
B07	0.44505530	0.30794600	Vitamin biosynthesis	0.13710930	Acclimation
B08	0.44618700	0.31972550	Aromatic compound biosynthesis	0.12646150	Acclimation
D09	0.25519710	0.13259720	Antibiotic degradation	0.12259990	Acclimation
B04	0.54481600	0.43492620	SCFA biosynthesis	0.10988980	Acclimation
D03	0.29173710	0.19845420	Sugar degradation	0.09328290	Acclimation
S02	0.19281000	0.15677840	Appendages	0.03603160	Acclimation
B10	0.05960097	0.03669565	Antibiotic biosynthesis	0.02290532	Acclimation

10.3.2 Warm population

10.3.2.1 Structural zeros

acclimation_samples <- sample_metadata %>%
  filter(Population == "Warm" & time_point == "Acclimation") %>% 
  dplyr::select(sample) %>%
  pull()

antibiotics_samples <- sample_metadata %>%
  filter(Population == "Warm" & time_point == "Antibiotics") %>% 
  dplyr::select(sample) %>% pull()

sample_sub <- sample_metadata %>%
  filter(Population == "Warm" & time_point %in% c("Acclimation", "Antibiotics"))

structural_zeros <- genome_counts_filt %>% 
  select(one_of(c("genome",sample_sub$sample))) %>%
  filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
  select(genome, where(~ is.numeric(.) && sum(.) > 0)) %>% 
   rowwise() %>% 
   mutate(all_zeros_acclimation = all(c_across(all_of(acclimation_samples)) == 0)) %>% 
   mutate(all_zeros_antibiotics = all(c_across(all_of(antibiotics_samples)) == 0)) %>% 
   mutate(average_acclimation = mean(c_across(all_of(acclimation_samples)), na.rm = TRUE)) %>% 
   mutate(average_antibiotics = mean(c_across(all_of(antibiotics_samples)), na.rm = TRUE)) %>% 
   filter(all_zeros_acclimation == TRUE || all_zeros_antibiotics==TRUE)  %>% 
   mutate(present = case_when(
      all_zeros_acclimation & !all_zeros_antibiotics ~ "antibiotics",
      !all_zeros_acclimation & all_zeros_antibiotics ~ "acclimation",
      !all_zeros_acclimation & !all_zeros_antibiotics ~ "None",
      TRUE ~ NA_character_
    )) %>%
   mutate(average = ifelse(present == "acclimation", average_acclimation, average_antibiotics)) %>%
   dplyr::select(genome, present, average) %>%
   left_join(genome_metadata, by=join_by(genome==genome)) %>%
   arrange(present,-average)

Acclimation

structural_zeros %>% 
  filter(present=="acclimation") %>% 
  count(phylum, name = "sample_count") %>%
  arrange(desc(sample_count)) 

# A tibble: 12 × 2
# Rowwise: 
   phylum               sample_count
   <chr>                       <int>
 1 p__Bacillota_A                 64
 2 p__Bacteroidota                32
 3 p__Bacillota                   16
 4 p__Cyanobacteriota              9
 5 p__Pseudomonadota               8
 6 p__Bacillota_B                  3
 7 p__Desulfobacterota             3
 8 p__Bacillota_C                  2
 9 p__Verrucomicrobiota            2
10 p__Actinomycetota               1
11 p__Campylobacterota             1
12 p__Elusimicrobiota              1

structural_zeros %>% 
  filter(present=="acclimation") %>% 
  select(1,5:10)

# A tibble: 142 × 7
# Rowwise: 
   genome                phylum               class               order                 family                 genus            species
   <chr>                 <chr>                <chr>               <chr>                 <chr>                  <chr>            <chr>  
 1 AH1_2nd_11:bin_000008 p__Verrucomicrobiota c__Verrucomicrobiae o__Verrucomicrobiales f__Akkermansiaceae     g__Akkermansia   s__    
 2 LI1_2nd_8:bin_000013  p__Verrucomicrobiota c__Verrucomicrobiae o__Verrucomicrobiales f__Akkermansiaceae     g__Akkermansia   s__    
 3 AH1_2nd_19:bin_000005 p__Bacillota_A       c__Clostridia       o__Lachnospirales     f__Lachnospiraceae     g__Roseburia     s__    
 4 AH1_2nd_20:bin_000009 p__Bacillota_A       c__Clostridia       o__Lachnospirales     f__Lachnospiraceae     g__Kineothrix    s__    
 5 AH1_2nd_18:bin_000011 p__Bacillota         c__Bacilli          o__Erysipelotrichales f__Coprobacillaceae    g__Coprobacillus s__    
 6 LI1_2nd_2:bin_000001  p__Bacillota_A       c__Clostridia       o__Oscillospirales    f__Ruminococcaceae     g__              s__    
 7 AH1_2nd_6:bin_000035  p__Bacillota         c__Bacilli          o__RF39               f__UBA660              g__Faecisoma     s__    
 8 LI1_2nd_4:bin_000003  p__Bacillota         c__Bacilli          o__Erysipelotrichales f__Erysipelotrichaceae g__Dielma        s__    
 9 AH1_2nd_1:bin_000039  p__Bacteroidota      c__Bacteroidia      o__Bacteroidales      f__Marinifilaceae      g__Odoribacter   s__    
10 AH1_2nd_2:bin_000018  p__Elusimicrobiota   c__Elusimicrobia    o__Elusimicrobiales   f__Elusimicrobiaceae   g__UBA1436       s__    
# ℹ 132 more rows

Antibiotics

structural_zeros %>% 
  filter(present=="antibiotics") %>% 
  select(1,5:10)

# A tibble: 7 × 7
# Rowwise: 
  genome                phylum            class                  order                 family                genus              species                       
  <chr>                 <chr>             <chr>                  <chr>                 <chr>                 <chr>              <chr>                         
1 AH1_2nd_7:bin_000003  p__Pseudomonadota c__Gammaproteobacteria o__Enterobacterales   f__Enterobacteriaceae g__Proteus         s__Proteus cibarius           
2 LI1_2nd_8:bin_000030  p__Actinomycetota c__Actinomycetia       o__Mycobacteriales    f__Mycobacteriaceae   g__Corynebacterium s__                           
3 LI1_2nd_8:bin_000077  p__Bacteroidota   c__Bacteroidia         o__Bacteroidales      f__Tannerellaceae     g__Parabacteroides s__Parabacteroides sp003480915
4 LI1_2nd_5:bin_000032  p__Bacillota_A    c__Clostridia          o__Christensenellales f__UBA1242            g__Caccovivens     s__                           
5 AH1_2nd_18:bin_000013 p__Bacteroidota   c__Bacteroidia         o__Bacteroidales      f__Tannerellaceae     g__Parabacteroides s__                           
6 LI1_2nd_5:bin_000069  p__Bacillota      c__Bacilli             o__RF39               f__UBA660             g__                s__                           
7 AH1_2nd_20:bin_000061 p__Bacillota      c__Bacilli             o__Lactobacillales    f__Enterococcaceae    g__Enterococcus    s__                           

Community elements differences in structural zeros

element_gift <- GIFTs_elements_community %>% 
  as.data.frame() %>% 
  rownames_to_column(., "sample") %>% 
  left_join(., sample_metadata[c(8,9)], by=join_by("sample"=="sample"))

uniqueGIFT_db<- unique(GIFT_db[c(2,4,5,6)]) %>% unite("Function",Function:Element, sep= "_", remove=FALSE)

element_gift %>%
    pivot_longer(-c(sample,time_point), names_to = "trait", values_to = "value") %>%
    group_by(trait) %>%
    summarise(p_value = wilcox.test(value ~ time_point)$p.value) %>%
    mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
    filter(p_adjust < 0.05)

# A tibble: 20 × 3
   trait p_value p_adjust
   <chr>   <dbl>    <dbl>
 1 B0220 0.00618   0.0497
 2 B0708 0.00280   0.0497
 3 B0709 0.00662   0.0497
 4 B1012 0.00618   0.0497
 5 D0203 0.00559   0.0497
 6 D0205 0.00280   0.0497
 7 D0206 0.00280   0.0497
 8 D0207 0.00280   0.0497
 9 D0210 0.00280   0.0497
10 D0211 0.00662   0.0497
11 D0213 0.00685   0.0497
12 D0305 0.00685   0.0497
13 D0308 0.00685   0.0497
14 D0610 0.00618   0.0497
15 D0706 0.00662   0.0497
16 D0902 0.00618   0.0497
17 D0905 0.00618   0.0497
18 D0908 0.00280   0.0497
19 D0911 0.00618   0.0497
20 S0104 0.00685   0.0497

10.3.2.2 Differential abundance analysis: composition

phylo_samples <- sample_metadata %>%
  filter(Population == "Warm" & time_point %in% c("Acclimation", "Antibiotics") )%>% 
  column_to_rownames("sample") %>% 
  sample_data()
phylo_genome <- genome_counts_filt %>% 
  filter(!genome %in% structural_zeros$genome) %>% 
  select(one_of(c("genome",rownames(phylo_samples)))) %>%
  filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
  select(genome, where(~ is.numeric(.) && sum(.) > 0)) %>%
  column_to_rownames("genome") %>% 
  otu_table(., taxa_are_rows = TRUE)
phylo_taxonomy <- genome_metadata %>% 
  filter(genome %in% rownames(phylo_genome)) %>% 
  column_to_rownames("genome") %>% 
  dplyr::select(domain,phylum,class,order,family,genus,species) %>% 
  as.matrix() %>% 
  tax_table()

physeq_genome_filtered <- phyloseq(phylo_genome, phylo_taxonomy, phylo_samples)

ancom_rand_output_2803 = ancombc2(data = physeq_genome_filtered, 
                  assay_name = "counts", 
                  tax_level = NULL,
                  fix_formula = "time_point",
                  rand_formula = "(1|individual)",
                  p_adj_method = "holm", 
                  pseudo_sens = TRUE,
                  prv_cut =0, 
                  lib_cut = 0, 
                  s0_perc = 0.05,
                  group = NULL, 
                  struc_zero = FALSE, 
                  neg_lb = FALSE,
                  alpha = 0.05, 
                  n_cl = 2, 
                  verbose = TRUE,
                  global = FALSE, 
                  pairwise = FALSE, 
                  dunnet = FALSE, 
                  trend = FALSE,
                  iter_control = list(tol = 1e-5, max_iter = 20, verbose = FALSE),
                  em_control = list(tol = 1e-5, max_iter = 100),
                  lme_control = lme4::lmerControl(),
                  mdfdr_control = list(fwer_ctrl_method = "holm", B = 100), 
                  trend_control = NULL)

save(ancom_rand_output_2803, 
     file = "data/ancombc_antib_accl_warm_2803.Rdata")

load("data/ancombc_antib_accl_warm_2803.Rdata")

taxonomy <- data.frame(physeq_genome_filtered@tax_table) %>%
  rownames_to_column(., "taxon") %>%
  mutate_at(vars(order, phylum, family, genus, species), ~ str_replace(., "[dpcofgs]__", ""))

ancombc_rand_table_mag <- ancom_rand_output_2803$res %>%
  dplyr::select(taxon, lfc_time_pointAntibiotics, p_time_pointAntibiotics) %>%
  filter(p_time_pointAntibiotics < 0.05) %>%
  dplyr::arrange(p_time_pointAntibiotics) %>%
  merge(., taxonomy, by="taxon") %>%
  mutate_at(vars(phylum, species), ~ str_replace(., "[dpcofgs]__", ""))%>%
  dplyr::arrange(lfc_time_pointAntibiotics)

colors_alphabetic <- read_tsv("https://raw.githubusercontent.com/earthhologenome/EHI_taxonomy_colour/main/ehi_phylum_colors.tsv") %>%
  mutate_at(vars(phylum), ~ str_replace(., "[dpcofgs]__", ""))  %>%
  right_join(taxonomy, by=join_by(phylum == phylum)) %>%
  dplyr::select(phylum, colors) %>%
  mutate(colors = str_c(colors, "80"))  %>% #add 80% alpha
    unique() %>%
    dplyr::arrange(phylum)

tax_table <- as.data.frame(unique(ancombc_rand_table_mag$phylum))
  
colnames(tax_table)[1] <- "phylum"
tax_color <- merge(tax_table, colors_alphabetic, by="phylum")%>%
    dplyr::arrange(phylum) %>%
    dplyr::select(colors) %>%
    pull()

ancombc_rand_table_mag%>%
      mutate(genome=factor(taxon,levels=ancombc_rand_table_mag$taxon)) %>%
ggplot(., aes(x=lfc_time_pointAntibiotics, y=forcats::fct_reorder(genome,lfc_time_pointAntibiotics), fill=phylum)) + 
  geom_col() + 
  scale_fill_manual(values=tax_color) + 
  geom_hline(yintercept=0) + 
  theme(panel.background = element_blank(),
        axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
        axis.text.x = element_text(size = 12),
        axis.text.y = element_text(size = 8),
        axis.title = element_text(size = 14, face = "bold"),
        legend.text = element_text(size = 12),
        legend.title = element_text(size = 14, face = "bold"),
        legend.position = "right", legend.box = "vertical")+
  xlab("log2FoldChange") + 
  ylab("Species")+
  guides(fill=guide_legend(title="Phylum"))

MAGs enriched in acclimation on the left side and in antibiotics on the right side

Phyla of the significant MAGs Acclimation

ancombc_rand_table_mag%>%
  filter(lfc_time_pointAntibiotics<0)  %>% 
  select(phylum, genus, species)

          phylum       genus                species
1   Bacteroidota Phocaeicola                       
2   Bacteroidota   Alistipes                       
3 Pseudomonadota  Salmonella    Salmonella enterica
4   Bacteroidota Bacteroides Bacteroides fragilis_B

Antibiotics

ancombc_rand_table_mag%>%
  filter(lfc_time_pointAntibiotics>0)  %>% 
  count(phylum, name = "sample_count") %>%
  arrange(desc(sample_count))  

          phylum sample_count
1   Bacteroidota            4
2      Bacillota            2
3    Bacillota_A            1
4    Chlamydiota            1
5 Pseudomonadota            1

ancombc_rand_table_mag%>%
  filter(lfc_time_pointAntibiotics>0)  %>% 
  select(phylum, genus, species)

          phylum           genus                  species
1   Bacteroidota Parabacteroides                         
2   Bacteroidota     Bacteroides                         
3   Bacteroidota     Bacteroides Bacteroides intestinalis
4    Bacillota_A      MGBC116941                         
5      Bacillota  Staphylococcus    Staphylococcus shinii
6      Bacillota      Ureaplasma                         
7   Bacteroidota     Bacteroides                         
8 Pseudomonadota                                         
9    Chlamydiota                                         

10.3.2.3 Differences in the functional capacity

sample_sub <- sample_metadata %>%
  filter(Population == "Warm" & time_point %in% c("Acclimation", "Antibiotics"))

genome_counts_filtered <- genome_counts_filt %>% 
    select(one_of(c("genome",sample_sub$sample)))

GIFTs_elements <- to.elements(genome_gifts,GIFT_db)
GIFTs_elements_filtered <- GIFTs_elements[rownames(GIFTs_elements) %in% genome_counts_filtered$genome,]
GIFTs_elements_filtered <- as.data.frame(GIFTs_elements_filtered) %>% 
  select_if(~ !is.numeric(.) || sum(.) != 0)

GIFTs_functions <- to.functions(GIFTs_elements_filtered,GIFT_db)

GIFTs_domains <- to.domains(GIFTs_functions,GIFT_db)

genome_counts_filtered_row <- genome_counts_filtered %>%
  mutate_at(vars(-genome),~./sum(.))  %>%
  column_to_rownames(., "genome") 

GIFTs_elements_community <- to.community(GIFTs_elements_filtered,genome_counts_filtered_row,GIFT_db)
GIFTs_functions_community <- to.community(GIFTs_functions,genome_counts_filtered_row,GIFT_db)
GIFTs_domains_community <- to.community(GIFTs_domains,genome_counts_filtered_row,GIFT_db)

10.3.2.3.1 Community elements differences:

MCI

GIFTs_elements_community %>%
  rowMeans() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = join_by(sample == sample)) %>%
  group_by(time_point) %>%
  summarise(MCI = mean(value), sd = sd(value))

# A tibble: 2 × 3
  time_point    MCI     sd
  <chr>       <dbl>  <dbl>
1 Acclimation 0.350 0.0225
2 Antibiotics 0.250 0.0665

MCI <- GIFTs_elements_community %>%
  rowMeans() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = join_by(sample == sample)) 

shapiro.test(MCI$value)

    Shapiro-Wilk normality test

data:  MCI$value
W = 0.88114, p-value = 0.03317

wilcox.test(value ~ time_point, data=MCI)

    Wilcoxon rank sum exact test

data:  value by time_point
W = 70, p-value = 0.0003291
alternative hypothesis: true location shift is not equal to 0

element_gift <- GIFTs_elements_community %>% 
  as.data.frame() %>% 
  rownames_to_column(., "sample") %>% 
  left_join(., sample_metadata[c(8,9)], by=join_by("sample"=="sample"))

uniqueGIFT_db<- unique(GIFT_db[c(2,4,5,6)]) %>% unite("Function",Function:Element, sep= "_", remove=FALSE)

significant_elements <- element_gift %>%
    pivot_longer(-c(sample,time_point), names_to = "trait", values_to = "value") %>%
    group_by(trait) %>%
    summarise(p_value = wilcox.test(value ~ time_point)$p.value) %>%
    mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
    filter(p_adjust < 0.05)%>%
  left_join(.,uniqueGIFT_db[c(1,3)],by = join_by(trait == Code_element))

element_gift_t <- element_gift  %>% 
  t() %>%
  row_to_names(row_number = 1) %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric)  %>%
  rownames_to_column(., "trait")

element_gift_filt <- subset(element_gift_t, trait %in% significant_elements$trait) %>% 
  t() %>%
  row_to_names(row_number = 1) %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric)  %>%
  rownames_to_column(., "sample")%>% 
  left_join(., sample_metadata[c(8,9)], by = join_by(sample == sample))

difference_table <- element_gift_filt %>%
  select(-sample) %>%
  group_by(time_point) %>%
  summarise(across(everything(), mean)) %>%
  t() %>%
  row_to_names(row_number = 1) %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric) %>%
  rownames_to_column(., "Elements") %>%
  left_join(.,uniqueGIFT_db[c(1,3,4)],by = join_by(Elements == Code_element)) %>% 
  arrange(Function) %>% 
  mutate(Difference=Acclimation-Antibiotics)%>% 
  mutate(group_color = ifelse(Difference <0, "Antibiotics","Acclimation")) 

difference_table %>%
  ggplot(aes(x=forcats::fct_reorder(Function,Difference), y=Difference, fill=group_color)) + 
  geom_col() +
  scale_fill_manual(values=c("#d57d2c","#003a45")) +
  geom_hline(yintercept=0) + 
  coord_flip()+
  theme(axis.text = element_text(size = 10),
        axis.title = element_text(size = 12),
        legend.position = "right", 
        panel.background = element_blank(),
          panel.grid.major = element_line(size = 0.15, linetype = 'solid',
                                colour = "grey"))+
  xlab("Microbial Functional Capacity") + 
  ylab("Mean difference")+
  labs(fill="Time point")

10.3.2.3.2 Community functions differences

MCI

GIFTs_functions_community %>%
  rowMeans() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = join_by(sample == sample)) %>%
  group_by(time_point) %>%
  summarise(MCI = mean(value), sd = sd(value))

# A tibble: 2 × 3
  time_point    MCI     sd
  <chr>       <dbl>  <dbl>
1 Acclimation 0.348 0.0158
2 Antibiotics 0.268 0.0585

MCI <- GIFTs_functions_community %>%
  rowMeans() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = join_by(sample == sample)) 

shapiro.test(MCI$value)

    Shapiro-Wilk normality test

data:  MCI$value
W = 0.85808, p-value = 0.01428

wilcox.test(value ~ time_point, data=MCI)

    Wilcoxon rank sum exact test

data:  value by time_point
W = 69, p-value = 0.0005759
alternative hypothesis: true location shift is not equal to 0

function_gift <- GIFTs_functions_community %>% 
  as.data.frame() %>% 
  rownames_to_column(., "sample") %>% 
  merge(., sample_metadata[c(8,9)], by="sample")

unique_funct_db<- GIFT_db[c(3,4,5)] %>% 
  distinct(Code_function, .keep_all = TRUE)

significant_functional <- function_gift %>%
  pivot_longer(-c(sample,time_point), names_to = "trait", values_to = "value") %>%
  group_by(trait) %>%
  summarise(p_value = wilcox.test(value ~ time_point)$p.value) %>%
  mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
  filter(p_adjust < 0.05)%>%
  left_join(.,unique_funct_db[c(1,3)],by = join_by(trait == Code_function))

Code_function	Acclimation	Antibiotics	Function	Difference	time_point
D02	0.42716670	0.21793910	Polysaccharide degradation	0.20922760	Acclimation
B02	0.62784110	0.47187490	Amino acid biosynthesis	0.15596620	Acclimation
D03	0.34000130	0.18813800	Sugar degradation	0.15186330	Acclimation
B08	0.44482920	0.33165890	Aromatic compound biosynthesis	0.11317030	Acclimation
D05	0.33300830	0.22491860	Amino acid degradation	0.10808970	Acclimation
D09	0.24738110	0.14240830	Antibiotic degradation	0.10497280	Acclimation
B06	0.70894950	0.61231610	Organic anion biosynthesis	0.09663340	Acclimation
B09	0.02947309	0.01568554	Metallophore biosynthesis	0.01378755	Acclimation